##############################
# STATE REACH AND DEVELOPMENT // PSRM // Replication
#
# Produces summary figure of robustness checks
#
##############################


# Load coefficients
all.coef.list <- readRDS("data/tempresults.rds")

# Get coefficient dataframe
coef.df <- do.call(rbind, lapply(names(all.coef.list), function(s){
  do.call(rbind, lapply(all.coef.list[[s]], function(m){
    data.frame(stub = s, 
               dep.var = m$dv,
               var = rownames(m$coef),
               coef = m$coef[,1],
               se = diag(m$clustervcv)^.5,
               stringsAsFactors = F)
  }))
}))

# Subset to main explanatory variables

## Adjust coef names for 1966
coef.df$var[coef.df$var =="log(1 + natcap.1966)" &  coef.df$stub == "1966roadnw"] <- "natcap.ln"
coef.df$var[coef.df$var =="log(1 + regcap.1966)" &  coef.df$stub == "1966roadnw"] <- "regcap.ln"

## Subset
coef.df <- coef.df[coef.df$var %in% c("natcap.ln", "regcap.ln"),]


# Sort plot
rob.spec <- list("tab1_main",
                 "taba7_migration",
                 "taba8_leads",
                 `Endogenous road changes` = c("roadscontr", "1966roadnw"),
                 "taba10_maaccess",
                 "taba11_ethnicpolitics",
                 `Capital controls` = "taba12_capitalcontr",
                 "popcontr")
plot.df <- do.call(rbind, lapply(1:length(rob.spec), function(i){
  df <- do.call(rbind, lapply(rob.spec[[i]], function(s){coef.df[coef.df$stub == s,]}))
  if(names(rob.spec)[i] != ""){
    df$type <- "plain"
    df <- smartbind(data.frame(stub = rep(paste0(names(rob.spec)[i],":"), length(unique(df$dep.var))),
                               dep.var = unique(df$dep.var), type = "bold",
                               stringsAsFactors = F),
                    df)
  } else {
    df$type <- "bold"
  }
  df
}))


# Rename Labels
plot.df$label <- plot.df$stub
plot.df$label[plot.df$label == "tab1_main"] <- "Baseline"
plot.df$label[plot.df$label == "taba10_maaccess"] <- "Market access controls"
plot.df$label[plot.df$label == "1966roadnw"] <- "stable 1966 road network" 
plot.df$label[plot.df$label == "taba11_ethnicpolitics"] <- "Ethnic politics controls" 
plot.df$label[plot.df$label == "roadscontr"] <- rep(c("add local road dummies", "add local road density"), each = 2)
plot.df$label[plot.df$label == "taba12_capitalcontr"] <- rep(c("add capital dummy", "add capital < 1hr dummy"), each = 2)
plot.df$label[plot.df$label == "popcontr"] <- "Local population controls"
plot.df$label[plot.df$label == "taba7_migration"] <- rep(c("Among non-migrants", "drop"), each = 2) ## drop second spec of these models, simple inteaction of controls with migrant dummy, does not give more insights. 
plot.df$label[plot.df$label == "taba8_leads"] <- "Controlling for pre-trends"
plot.df <- plot.df[plot.df$label != "drop",]

# ... to factors and numeric positions
plot.df$label <- factor(plot.df$label, levels = rev(unique(plot.df$label)), ordered = T)
plot.df$pos <- as.numeric(plot.df$label)

# Rename variables
all.dv <- unlist(lapply(main.data.names, get_depvar))
names(all.dv) <- gsub("(","\n(", names(all.dv), fixed = T)
plot.df$outcome.name <- factor(names(all.dv)[unlist(lapply(plot.df$dep.var, function(d){which(all.dv == d)}))], 
                               levels = names(all.dv), ordered = T)
plot.df$var.name <- ifelse(grepl("natcap", plot.df$var), "National cap.", "Regional cap.")
plot.df$var.name <- factor(plot.df$var.name, levels = unique(plot.df$var.name), ordered = T)


# ... print
png(file.path(fig.path, paste0("fig5_robchecks.png")), width = 10, height = 5, units = "in", res = 300)
p <- ggplot(plot.df, aes(x = coef, y = pos)) +
  geom_vline(xintercept = 0, lty = 2, col = "grey") +
  geom_point() +
  geom_segment(aes(x = coef + se * 1.96,
                   xend = coef - se * 1.96,
                   y = pos, yend = pos)) +
  facet_wrap(~ outcome.name + var.name, nrow = 1, scales = "free_x") +
  xlab("Marginal effect of time to national/regional capital (log)") + 
  ylab(paste0("")) + 
  scale_y_continuous(breaks = unique(plot.df$pos), labels = unique(plot.df$label)) +
  scale_x_continuous(breaks = scales::pretty_breaks(n = 3)) +
  theme_minimal()  +
  theme(axis.text.y = element_text(face = unique(plot.df[,c("label","type")])$type, 
                                   size = 10),
        strip.text = element_text(size=8),
        panel.spacing = unit(1, "lines")) +
  NULL
print(p)
dev.off()

